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We apply the consistent discretization scheme to general relativity particularized to the Gowdy 
space-times. This is the first time the framework has been applied in detail in a non-linear generally- 
covariant gravitational situation with local degrees of freedom. We show that the scheme can be 
correctly used to numerically evolve the space-times. We show that the resulting numerical schemes 
are convergent and preserve approximately the constraints as expected. 



--- I. INTRODUCTION 
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, The "consistent discretization" approach [Ij has proven quite attractive to deal with conceptual issues in quantum 

' gravity P|. In this approach one approximates general relativity with a discrete theory that is constraint-free. As a 
consequence, it is free of the hard conceptual issues that plague quantum gravity. One can, for instance, solve the 

' "problem of time" and one discovers that there is a fundamental decoherence 0, Q induced in quantum states that 
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can yield the black hole information puzzle unobservable 

A looming question in this approach is how well can the constraint-free discrete theory approximate general rela- 
tivity? In this approach one starts with a Lagrangian, discretizes space-time and then works out discrete equations of 
motion that determine not just the usual variables of general relativity but the lapse and the shift as well. This has 
the attractive aspect of yielding a consistent set of discrete equations (they can all be solved simultaneously, unlike in 
usual discretizations where the constraints are not preserved upon evolution; it should be emphasized that preserving 
the constraints is becoming a central focus in state-of-the-art numerical relativity.) However, the resulting equations 
for determining the lapse and the shift are non-polynomial and of a rather high degree. Since these equations have 
no counterpart in the continuum theory, one does not have an underlying mathematical theory to analyze them as 
one has in the case of the other discrete equations. For instance, it is not obvious that the solutions of the equa- 
. tions will be real. Or that they do not oscillate wildly upon each step of evolution. We had probed these issues in 
some detail in some cosmological models There the consistent discretization approach works well, approximating 
the continuum theory in a controlled fashion and covering all of the phase space. But in the cosmological case the 
O I equations also simplify significantly. Moreover, the diffeomorphism constraint does not play a role. In addition to 
this, it is known that discretizations of evolution equations can pose unique challenges. For instance, it is not clear 
5h . that the resulting scheme produced by the consistent discretizations is hyperbolic in any given sense. With regular 
^f^' discretizations, the use of formulations that are not hyperbolic has led to problems. The type of discretizations that 
^ I arise in the consistent discretization approach are also "forward in time centered in space" . This is due to the fact 
' that the use of centered derivatives in time complicates the construction of the canonical formulation that is central 
r\ to the use of consistent discretizations in quantum gravity. All this raised significant suspiciousness about how well 
' these discretizations could approximate the continuum theory. An encouraging point was that at least for linearized 
gravity, the consistent discretizations yield a "mimetic" discretization ^ that appears to be stable. But the linearized 
case has many peculiarities that do not carry over to the non-linear case (namely the mimetism). It is also true 
that at the classical level the "consistent discretization" approach may be considered as an extension of the idea of 
variational integrators to singular Lagrangians (systems with constraints treated in the Dirac fashion). Variational 
integrators have proven to have many desirable properties, at least for unconstrained systems (see |lOj for a recent 
review), but they have only been applied in a limited way to systems with constraints (only holonomic or at best 
constraints linear in the momenta [ll|) without treating them in the usual Dirac fashion. 

It is clear that the viability of the approach has to be probed in situations with field-theoretical degrees of freedom. 
A situation of this type arises in the Gowdy cosmologies J^J . These are spatially compact space-times with two space- 
like commuting Killing vector fields. This allows the introduction of coordinates where the metric depends only on time 
and an angular variable which we will call 6. The consistent discretization approach applied to this model yields quite 
involved equations, as we will see. They cannot be solved analytically and have to be tackled numerically. Moreover, 
since we are dealing with non-linear coupled algebraic systems, the solution has to be constructed approximately, and 
the resulting systems are rather large. In addition to this, the presence of global constraints derived from the compact 
topology yields some of the systems unsolvable directly and this issue has to be addressed again by an approximate 
technique. 

It should be made clear from the outset that we are not at all interested in being competitive in the well developed 
field of numerical simulations of Gowdy cosmologies (see Berger for a review). Usually numerical simulations of 
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Gowdy take advantages of special gauges and other features that we will not take into account in our approach. 

To simplify matters we will concentrate on the "single polarization" Gowdy case and we will further restrict 
to a particular subcase that can be seen as a "minisuperspace" by setting one of the variables and its conjugate 
momentum to zero. From the point of view of the evolution equations and the dynamics in terms of coordinate 
dependent quantities, the complexity of the equations is virtually similar to the full single-polarized case, one just 
is faced with a smaller number of equations. Moreover, if one were to write the full one polarization case and take 
suitably restricted initial data, the evolution is preserved exactly within the minisuperspace considered. 

The organization of this article is as follows. In the next section we briefly review the formulation of Gowdy 
cosmologies mostly to fix notation. We then proceed to work out the consistent discretization. In section III we 
discuss the numerical scheme and present the numerical results. 

II. GOWDY COSMOLOGIES 

A. Continuum formulation 

Gowdy cosmologies are space-times with two space-like commuting Killing vector fields. We take the topology of 
the spatial slices to be a three-torus T^. Following Misner we parameterize the spatial metric by, 

ds^ = e-^-^^/'^d9' + e'^ {e^da^ + e~^dd') (1) 

where functions r. A, (3 depend on time t and on one of the spatial variables which we call 9 e [0, 27r]. The other two 
spatial variables are a and S and will play no role from now on. The Hilbert action particularized for these kinds of 
metrics reads, 

S^^JdtJde^pxX+PrT+ppP-N^C'^,) (2) 

where N^, /i = 0, 1 are related to the lapse and the single remaining component of the shift (Lagrange multipliers) by 
exp(— 3r/2 -I- A/4), Ng = Ng exp (t -|- A/2). Also are related to the Hamiltonian and the single 
remaining component of the momentum constraint by the obvious rescalings, 

C' = 4p'^+PrT' +ppf3' +pxX'. (4) 

We have also chosen the irrelevant spatial coordinates in such a way that J da J dS — 8. 

The theory can be consistently reduced by setting P = pp — 0. This corresponds to a one-parameter family of 
space-times. Since in our approach the coordinates are not fixed, the treatment is highly non-trivial, and the problem 
appears to be "infinite dimensional" (although one is strictly speaking dealing with a zero dimensional situation, as 
is the case in minisuperspaces). The exact solution of the evolution equations for this case can be written, in a given 
gauge choice (as shown by Misner). If one chooses lapse equal one, shift equal zero, r' = 0, (px)' — one completely 
fixes gauge and the solution results r = ct, p\ = c, Pr = and A = 0. We will see that in our approach we can choose 
initial data where the shift approximately vanishes upon evolution. This is good since it will allow us to compare 
scalars computed in the numerical solution with those of the exact solution. If the shift is non-vanishing then the 
situation complicates since one does not know at which points to compare the invariants (one is facing the classic 
"metric equivalence problem" , i.e., comparing metrics in two different coordinate systems to see if they are the same). 

B. Consistent discretization 

We now proceed to work out the consistent discretization of the theory. We assume that the time coordinate is 
related to a discrete variable n hy t = nAt, n = . . . N. The spatial coordinate 9 is also discretized 9 — ■mA9 with 
m = . . . M. The action then becomes, 

S = ^i(n,n + l) 

N M ^ 

— p\{n,m) {X{n + l,m) — X{n,m)) + pr{n,rn) (T{n + l,Tn) — T{n,m)) (5) 

n—Q m— 
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-M{n, m) 



Px{n, m)pr{n, m) + e'*'^^"'™) ^4T(n, m + 1) + 4T(r7,, m - 1) - 8T(n, m) 



+8 (T(ri, m + 1) — T(n, m — 1))^ + (T(n, to + 1) — T(n, to — 1)) (A(n, to + 1) — A(n, m — 1)) 



-N{n, m) 



4:p\{n, m + 1) — 4pa('t-, to — 1) 



+p\{n, to) (A(n, TO + 1) — A(ri, m — 1)) + Pt('t-, ™) ('''(f^j m + 1) — T(n, to — 1)) 

where we have rescaled the momenta p(n,m) = p{t,9)A9, and rescaled and relabeled the shift N{n,m) = 
Ng{t,e)At/{Ae), and the lapse M{n,m) = N{t,e)At/{A0). The reader should keep in mind that in the rest of 
this paper when we refer to the "lapse" we are really referring at the lapse times the time interval in the discretization 
divided by the spatial interval. For example, we will therefore encounter statements like "making the lapse small" 
signifying that the discretization step, measured in an invariant fashion, is becoming smaller. 

We now proceed to define the canonical variables from the action. The reader should not be confused by the 
use of variables named "p" in the action; in the consistent discretization scheme all variables are initially treated as 
configuration variables (i.e. one is working in a first order formulation of the theory, see for a full discussion.) 
One defines canonical momenta at instant n and at instant n + I. Let us consider for instance the variable A. We 
define the canonical momenta at instant n + 1 by, 

dL 

+ ^' " 5A(n + l,TO) = 
and at instant n + I, using the Lagrange equations of motion one has. 



dL 



dX{n, to) 



P\{n, m) + M{n, m — 1) (T(n, m) — T(n, to — 2)) — M{n, to + 1) (r^n, to + 2) — T(n, to)) 
+N{n, m —- l)px{n, m — 1) + N{n, m + \-)p\{n^ to + 1). (7) 



One can combine these equations to eliminate the variable px{n, ni) from the model and yield an evolution equation 
in terms of variables that are genuinely canonically conjugate, 

P^{n + 1, to) = m) - A/(n, m)e'^^("'™' (T(n, to) - T(n, to + 1)) - M(n, m - l)e4^(«^™-i) (r(n, m) - r(n, m - 1)) 

+N{n, m)P^{n + 1, to) - N{n, to - l)P^{n + 1, to - 1). (8) 

A similar procedure yields an (implicit) evolution equation for P'^, 

P^{n+l,m) = P''{n,m) 

—M{n, TO ) (^4e4^("'™) (irin, to + 1) + 4T(n, to - 1) - 8T(n, to) + 8 (T(n, to + 1) - T(n, m) f 

+ {T{n, m + 1) — T{n, m)) (A(n, m + 1) — A(n, to)) ^ 

_^4g4r(n,m) (-_g _ j^g^^^^ TO + 1) + 16T(n, to) + A(n, to) — A(n, to + 1)) ) 

+N{n, m)P^{n + 1, to) - N{n, m - l)P^{n + 1, to - 1) - 4M(n, to + i)e4^(«.™+i) 

-Af(n, TO - i)e4^("'™"i) (4 + 16r(n, m) - 16T(n, m - 1) + A(n, m) - A(n, m - 1)) , (9) 

and for A, 

\{n + 1,to) = X{n,m) + M{n,m)P''{n + 1,to) + N{n,m) (-4 + X{n,m+ 1) - X{n,m)) + 4N{n,m- 1) (10) 
and T, 

t(7t, + 1, to) = r(n, to) + AI{n, m)P'^(n + 1, to) + N(n, m) {T{n, to + 1) — T{n, m)) . (11) 

A similar treatment for the Lagrange multipliers (defining their canonical momenta and combining the equations 
at n and n + 1) yields the "pseudoconstraints" , 

4P^{n + 1, TO + 1) - AP^in + 1, to) 

+P^{n + 1, to) (A(n, TO + 1) — A(n, m)) + P^ (n + 1, to) (T(n, to + 1) — riji, to)) = (12) 
P^{n + 1, m)P^{n + 1, to) + e''^("^™) [4T(n, m + 1) + 4T(n, m - 1) - 8r(n, to) 



+8 (T(n, m + 1) — T(n, m))^ + (T(n, to + 1) — T(n, m)) (A(n, to + 1) — A(n, m)) 



0. 



(13) 
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The terminology "pseudoconstraints" reflects the fact that these equations can be seen as discretizations of the 
constraint equations of the continuum theory. However, they are not constraints for the discrete theory in the usual 
sense of the word since they involve variables at different instants of time. 

The set of equations H8I13|I constitute six nonlinear coupled algebraic equations for the six unknowns, 
T, X, P'^ , P"^, N, M, that is, they determine the canonical variables and the Lagrange multipliers. In fact, this is 
an oversimplified view of the situation, since the equations link variables at different spatial points m. In order to 
solve them, one has to take into account that all variables in the Gowdy problem are periodic, for instance in a lattice 
with mm points, r(0) = T(mm), etc. and then one is left with a system of 6 x mm equations with 6 x mm unknowns. 
So if one decides to use 10 spatial points these are 60 equations with 60 unknowns. 

In fact, one does not need to tackle the system in an entirely coupled fashion, since the two evolution equations for 
A and t are explicit. We will discuss the solution strategy in more detail in the next section. 

A final remark is that Van Putten proposed a scheme which also can be characterized as determining the lapse and 
shift, but it differs significant from ours (in particular it is not motivated variationally) . He also applied it to the 
Gowdy cosmologies [l6.j . 

III. NUMERICAL RESULTS 
A. Choice of initial data 

In order to evolve the equations we need to set initial data. The choice of initial data is somewhat delicate, since it 
will determine not only the solution but also the lapse and the shift. Although one can give arbitrary initial data, it 
would be desirable if the resulting lapse were of the same sign across the spatial manifold. If one does not make this 
choice, the scheme is able to evolve, but portions of the spacetime will be evolved forward in time and some portions 
it will evolve backward in time. There is nothing wrong with this, but it is not what is traditionally considered in 
numerical evolutions. We may also wish to choose a vanishing shift, although this is not important. 

In order to ensure that the lapse and shift are approximately in line with what we desire, we proceed in the following 
way. We start with some initial data for the variables t(0), A(0), iV(0), M(0). We use equations 112I13|) to determine 
P'^(l), P'^(l) and then we use equations (|8I9|I to determine P'^(O), P^(0). This completes the initial data set. That 
is, we have given us r, A and the lapse and the shift that we desire, and this determined the canonically conjugate 
momenta. 

At this point the reader may be confused. Usually in numerical relativity one gives the metric and its canonical 
momenta at a given hypersurface, and they are chosen in such a way that they satisfy the constraints. To implement 
the consistent discretization evolution in this case one would consider the metric as evaluated at and the canonical 
momenta evaluated at 1. One then chooses lapse and shift at and evolves. Presumably one would choose a small 
lapse in order to generate an evolution that is close to the continuum theory (recall that we refer here to the rescaled 
lapse). One could do things the other way around, choosing the metric at 1 and the canonical momentum at and 
evolve backwards. The choice we made to of initial data is equivalent to finding a solution of the usual constraints at 
a given hypersurface with A, r given. 

B. Evolution 

We need to give a prescription that, given T(n), A(n), P'^(n), P^{n) determines the variables at instance n + 1 and 
in the process determines the Lagrange multipliers. One could take the complete set of equations and solve them 
simultaneously. We did not proceed like this. Since this was the first exploration of the system we thought it would be 
more instructive and offer more control on the situation to solve the equations in sub-systems. We will see however, 
that in the end one pays a price for this. 

The scheme we used to solve the system is as follows: using equations H12I13|) we determine P'^(n-t-l), P'^(n-l-l). We 
then take these values together with A(n), r(n), P'^(n), P'^(n) and use the equations (|8I9|I to determine the Lagrange 
multipliers N{n),M{n). With these multipliers and the initial values, it is possible, using H10lll|) to determine 
A(n + 1), r(n + 1), which would complete the evolution process from level n to n + 1. 

Unfortunately, there is a problem with this procedure. When we try to use equations H8I9|) to determine M{n), N{n) 
one finds that the system is indeterminate. One can determine all the Lagrange multipliers except one, which we 
can choose to be Ml = M{n,m = 1). This has to do with a symmetry of the Gowdy problem. The action (0) only 
depends on derivatives of the variable A, therefore it is invariant under the addition of a constant to that variable. 
Using Noether's theorem one can find out that there is a conserved quantity, J P^{9)d9. A similar result holds in the 

discrete action, the conserved quantity becomes X]m=o^^('^)- '-"^'^ imposes that this quantity be conserved, this 
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can be used to determine entirely the Lagrange multipliers. Another way of seeing this is to notice that one can cast 
the conservation law as a constraint on the initial data, 



m=0 



M 

^ P^{72 + l,m)[\{n,m),T{n,m)] = ^ P^{n,m) 

m=0 



(14) 



that is, writing the P'*' at n + 1 as a function of the initial data and imposing the conservation one has a constraint 
on the initial data. Imposition of this constraint is enough to determine entirely the Lagrange multipliers. 

This is a pathology of our choice of scheme of solution of the coupled system of equations, dealing with it by parts. 
If we had chosen to solve the entire non-linear system H8I11|I such system would have been well defined and would 
have automatically preserved the conserved quantity upon evolution. This would have, however, forced us to deal 
with a significantly larger coupled non-linear system. 

Numerically, the way we will implement the conservation is to choose arbitrarily the value of one component of the 
lapse. We then run the evolution scheme we outlined and we will check if X]m=o varied. If it has, we will go 

back and adjust the arbitrary component of the lapse until the conserved quantity is kept constant. We will do this 
via a Newton-Raphson technique. The fiow diagram of the logic is shown in figure 1. For the solution of the various 
non-linear systems we discussed we use the routine TENSOLVE (T7|. 



n 



n+1 (10,11-)-'^ ^ 



n+2 



(12,13) 



±1 



N(M1) M(M1) 



Ml 



M 



M 



Conserved? 



FIG. 1: The logic of the evolution scheme. 



C. Results for the dynamical variables 

In figure|21we show the variable A as a function of the spatial points and as a function of time, for a simulation with 
8 and 40 spatial points. The time is measured by the variable r, which is a good measure since it is invariant under 
coordinate changes of the t,6 variables. Since we have chosen the shift to be close to zero, and it is preserved that 
way upon evolution as we shall see, the value of the angle 6 can be taken to be an invariant. We see that the variable 
A exhibits the types of oscillations of increasing frequency that one encounters in the exact solution for the Gowdy 
spacetime. The evolutions with different spatial resolutions all have the same initial data for A(0), t(0), Af(0), iV(0). 
The explicit form of initial data chosen are, 

T(0,i) = -0.5 + 0.01 sin(27ri/M), (15) 

A(0,i) = 0.001 sin(27ri/M) -I- 0.0025 sin(47ri/M), (16) 

M{0,i) = 0.0005, (17) 

N{0,i) = 10"'^sin(27ri/Af). (18) 

Following the initial data construction we outlined before, this produces differing values of P^{0), P'^(O) depending 
on the number of points of the spatial grid. As shown in figure El there is convergence in the form of the resulting 
P'''(0), P'^(O) when one refines the spatial grid. 

The A(0), t(0), P'*'(0), P'^(O) produced by the initial data solving procedure, in the limit in which one makes infinitely 
large the number of spatial points, is not a solution of the initial value constraints of the continuum theory. However 
one can approximate a solution of the constraints of the continuum theory arbitrarily by choosing the value of the 
initial (rescaled) lapse to be very small. 
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40 points 
S points 




FIG. 2: The variable A as a function of space (6) and time, measured by the variable r. We show two resolutions, one with 
8 and another with 40 spatial points. The latter represents what is practical to run on a workstation today. One sees that 
the two resolutions track each other well for a while, but things deteriorate as time increases (see for instance the third line of 
crosses from the right). The run with 8 spatial points just does not have enough resolution to track the features of the solution 
in question, as can be seen from the figure. 




FIG. 3: The variable in the initial data converges to a continuum solution of the pseudo-constraints, and would converge 
to a solution of the continuum usual constraints of GR if the choice of initial data for the lapse and shift were zero. 



Figure^ shows the value of the average of the shift across the grid as a function of time. We see that in aU cases 
the shift is small, and is smaller for higher resolutions. This is desirable in order to compare with the exact solution 
which is in a gauge with zero shift. 

Figure [S] shows the Riemann tensor for three resolutions. More precisely, it represents component i?3434 (where 
3, 4 represent the ignorable coordinates) as a function of r for 9 = tt. The reason we chose this component is that 
it behaves as a scalar with respect to t,6 diffeomorphisms. Since we determine the lapse and the shift dynamically, 
for different resolutions we have different lapses and shifts and therefore different coordinates. As a consequence 
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FIG. 4: The L2 norm of the shift, which shows we are approximately in a gauge with zero shift. 

one cannot easily compare the other components of the Riemann tensor for different resolutions. One sees that the 
method converges, i.e. for smaller spatial separations the value of the Riemann tensor approximates well that of flat 
space-time. 




FIG. 5: Convergence of the Riemann tensor for three different resolutions. Since we are working with an almost vanishing shift, 
the components of the Riemann tensor can be treated as observables and compared among different resolutions directly. With 
a non- vanishing dynamically determined shift this would not be possible. In such cases one can only compare observables of 
the theory. 

We know the exact form of the Riemann tensor. In figureElwe use it to evaluate the relative error in the evaluation. 
The exact form is i?3434 = c^e""^ where c is the only parameter present in the minisuperspace (recall t = ct). To 
avoid having to determine the constant we actually plot in figure the following quantity, 

fi3434(T) _ g-(r-r(l)) 
-R3434(t(1)) (^Q\ 

which yields the relative error. 

In figure [3 we show the L2 norm of the square root of the sum of the squares of the constraints of the continuum 
theory evaluated in the discrete theory, for three resolutions. As expected, the magnitude of the constraints correlates 
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FIG. 6: Relative error of the Riemann tensor compared to the exact solution. 
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FIG. 7: L2 norm of the constraints of the continuum theory evaluated in the discrete theory. As can be seen they are preserved 
well and in a convergent fashion. As argued in the text, the value of the constraints is a measure of the error of the evolutions 
in the discrete theory and this plot confirms it. 



well with the error in the evolution scheme. We are displaying the values of the constraints with any type of 
normalization. This may be misleading. As we see, although there is convergence there is also exponential growth. 
The growth can be adscribed to the general exponential growth of variables in Gowdy, particularly the factor exp(4T) 
that appears in the Hamiltonian constraint. To compensate for this in figure |S1 we display the same data divided by 
the offending factor, and one cannot distinguish constraint growth. 



D. The flat sector 



If one chooses initial data for the evolution equations 1)819110111(1 such that — and r = 0, one can check 
that, in the continuum theory, one is producing flat space-time. We would like to check if the discrete theory is 
able to reproduce, at least in an approximate way, this behavior. This is usually known as the "gauge wave test" as 
is described, for instance in the "apples with apples" (www.appleswithapples.org) project, and it is known to be a 
somewhat significant hurdle for codes to pass (at least in full 3D). 

In the discrete theory the pseudo-constraints (|12I13|) are identically satisfied and the lapse and shift are free. Let 
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FIG. 8: L2 norm of the constraints of the continuum theory evaluated in the discrete theory normahzed by the factor exp(4r) 
that appears in the Hamiltonian constraint. One sees convergence with absence of growth in the normahzed constraints. The 
reader may be puzzled as to why the curves do not start at the same point. The reason is that we are displaying the curves 
starting at the first iteration, not at the initial data. For very low resolutions the variables change quite a bit in the first 
iteration. 



us start by considering a slicing with M = 1, iV = 0. There are only two evolution equations left, 

P^{n + l,m) = P^(n,m) + A(n,m + 1) - 2A(n,TO) + A(n,m - 1) (20) 
A(n + l,m) = P^(n + l,m) + A(n,m). (21) 

These equations can be combined to yield a single equation for A, 

A(n + 1, m) = A(n, rn + 1) + A(n, m - 1) - A(n - 1, to), (22) 

and it is remarkable to notice that the exact solution of the discrete equations for A is given by a plane wave, 

A(n, to) — f{n + to) + g{n — to), (23) 

and therefore the space-time metric is manifestly a "gauge wave" . Also remarkable is that if one computes the discrete 
Riemann tensor all its components vanish identically. In this sense, our formalism aces one of the "apples with apples" 
tests exactly. If we choose M = constant instead of one we will get boosted slices and again we will get the gauge 
wave as an exact solution and a vanishing Riemann tensor. If we choose more complicated slices, then the solution 
will be reproduced approximately. The resulting equations are, 

P^(n + 1,to) - P''(n, to) = A/(n, to)(-8 - \{n,m + 1) + A(n,m)) - M{n,m - + X{n,m) - \{n,ni - 1)) 

+N{n, 'm)P^{n + 1, to) - N{n, m - l)F^(n + 1, to - 1) - 4M (n, to + 1) (24) 
A(n + 1, to) - A(n, to) = A/(n, TO)P^(n + 1, to) + A^(n, to)(-4 + A(n, to + 1) - A(n, to) + 4A^(n, to - 1) (25) 

We have studied this case numerically and the code proves convergent and long term stable (we could not find 
evolutions that crashed) provided that one limits oneself to values of M less than unity, otherwise one can explicitly 
check that the evolution system has eigenvalues of the amplification matrix larger than one (this is just the Courant 
condition). 



IV. DISCUSSION 



We have applied for the first time the "consistent discretization" approach to a nonlinear situation with field 
theoretic degrees of freedom, the Gowdy one polarization cosmologies. We see that one can successfully evolve in a 
convergent and stable fashion for about one crossing time. All of the expected features of the approach are present. 
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One encounters that the Lagrange multipliers are determined by the evolution and that at some points they can 
become complex. Evolution can be continued by backtracking and switching to a different root of the equations that 
determine the Lagrange multipliers. We have not resorted to this type of technique for the results presented up to now 
in this paper. However, we have been able to extend the runs to about ten crossing times (r = 2.8) using these types 
of techniques, but we have not carried out detailed convergence studies. The codes finally stop due to exponential 
overflows in the resolution of the linear system derived from the implicit nature of the equations. 

The evaluation of the constraints of the continuum theory indicate that they are well preserved in a convergent 
fashion and can be used as independent error estimates of the solutions produced. 

It should be noted that one has choices when one discretizes the theory. We have chosen here to discretize the full 
model without any gauge fixing. One could have, however, resorted to fixing gauge and discretizing a posteriori. There 
are various possibilities, since one can also partially gauge fix and discretize. Partial gauge fixings may be desirable 
in realistic numerical simulations in order to have some control on the gauge in which the simulations are performed. 
For instance, in binary black hole situations it has been noted that the use of corrotating shifts is desirable. One 
could envision a scheme in which one gauge fixes the shift to the desired one and uses the consistent discretization to 
determine the lapse. 

The Gowdy example unfortunately has limitations as a test of the scheme given the exponential nature of the metric 
components. This prevents any scheme from achieving long term evolutions. A scheme like ours that involves complex 
operations involving many powers of the metric components is even more limited by this problem. We have tried to 
improve things by using quadruple precision in our Fortran code, but since one is battling exponential growth, this 
only offers limited help. 

Summarizing, the scheme works, but given the limitations of the example it is not yet clear that it will prove useful 
in long term evolutions in other contexts. Turning to quantum issues, the complexity exhibited in the determination 
of the lapse and shift suggests that quantization of these models will have to be tackled numerically. This paper can 
be seen as a first step in this direction as well. 
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